home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / README < prev    next >
Text File  |  1988-07-11  |  13KB  |  268 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that this notice is preserved and that due credit is given
  7.  * to the University of California at Berkeley. The name of the University
  8.  * may not be used to endorse or promote products derived from this
  9.  * software without specific prior written permission. This software
  10.  * is provided ``as is'' without express or implied warranty.
  11.  *
  12.  * All recipients should regard themselves as participants in an ongoing
  13.  * research project and hence should feel obligated to report their
  14.  * experiences (good or bad) with these elementary function codes, using
  15.  * the sendbug(8) program, to the authors.
  16.  *
  17.  * K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan.
  18.  * Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85.
  19.  *
  20.  *    @(#)README    5.2 (Berkeley) 4/29/88
  21.  */
  22.  
  23. ******************************************************************************
  24. *  This is a description of the upgraded elementary functions (listed in 1). *
  25. *  Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over    *
  26. *  from 4.2BSD without change except perhaps for the way floating point      *
  27. *  exception is signaled on a VAX.  Three lines that contain "errno" in erf.c*
  28. *  (error functions erf, erfc) have been deleted to prevent overriding the   *
  29. *  system "errno".                                                           *
  30. ******************************************************************************
  31.  
  32.  
  33. 0. Total number of files: 40 (Note for Sprite users:  the distribution of
  34.     files among directories has been changed slightly from what is
  35.     listed below in order to conform to Sprite's standard directory
  36.     organization).
  37.  
  38.         IEEE/Makefile   VAX/Makefile    VAX/support.s   erf.c       lgamma.c
  39.         IEEE/atan2.c    VAX/argred.s    VAX/tan.s       exp.c       log.c
  40.         IEEE/cabs.c     VAX/atan2.s     acosh.c         exp__E.c    log10.c
  41.         IEEE/cbrt.c     VAX/cabs.s      asincos.c       expm1.c     log1p.c
  42.         IEEE/support.c  VAX/cbrt.s      asinh.c         floor.c     log__L.c
  43.         IEEE/trig.c     VAX/infnan.s    atan.c          j0.c        pow.c
  44.         Makefile        VAX/sincos.s    atanh.c         j1.c        sinh.c
  45.         README          VAX/sqrt.s      cosh.c          jn.c        tanh.c
  46.  
  47. 1. Functions implemented :
  48.     (A). Standard elementary functions (total 22) :
  49.         acos(x)                 ...in file  asincos.c 
  50.         asin(x)                 ...in file  asincos.c
  51.         atan(x)                 ...in file  atan.c
  52.         atan2(x,y)              ...in files IEEE/atan2.c, VAX/atan2.s
  53.         sin(x)                  ...in files IEEE/trig.c,  VAX/sincos.s
  54.         cos(x)                  ...in files IEEE/trig.c,  VAX/sincos.s
  55.         tan(x)                  ...in files IEEE/trig.c,  VAX/tan.s
  56.         cabs(x,y)               ...in files IEEE/cabs.c,  VAX/cabs.s
  57.         hypot(x,y)              ...in files IEEE/cabs.c,  VAX/cabs.s
  58.         cbrt(x)                 ...in files IEEE/cbrt.c,  VAX/cbrt.s
  59.         exp(x)                  ...in file  exp.c
  60.         expm1(x):=exp(x)-1      ...in file  expm1.c
  61.         log(x)                  ...in file  log.c
  62.         log10(x)                ...in file  log10.c
  63.         log1p(x):=log(1+x)      ...in file  log1p.c
  64.         pow(x,y)                ...in file  pow.c
  65.         sinh(x)                 ...in file  sinh.c
  66.         cosh(x)                 ...in file  cosh.c
  67.         tanh(x)                 ...in file  tanh.c
  68.         asinh(x)                ...in file  asinh.c
  69.         acosh(x)                ...in file  acosh.c
  70.         atanh(x)                ...in file  atanh.c
  71.                         
  72.     (B). Kernel functions :
  73.         exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh
  74.         log__L(s)   ...in file log__L.c, used by log1p/log/pow
  75.         libm$argred ...in file VAX/argred.s, used by VAX version of sin/cos/tan
  76.  
  77.     (C). System supported functions :
  78.         sqrt()      ...in files IEEE/support.c, VAX/sqrt.s
  79.         drem()      ...in files IEEE/support.c, VAX/support.s
  80.         finite()    ...in files IEEE/support.c, VAX/support.s
  81.         logb()      ...in files IEEE/support.c, VAX/support.s
  82.         scalb()     ...in files IEEE/support.c, VAX/support.s
  83.         copysign()  ...in files IEEE/support.c, VAX/support.s
  84.         rint()      ...in file  floor.c
  85.  
  86.  
  87.    Notes: 
  88.        i. The codes in files ending with ".s" are written in VAX assembly 
  89.           language. They are intended for VAX computers.
  90.  
  91.           Files that end with ".c" are written in C. They are intended
  92.           for either a VAX or a machine that conforms to the IEEE 
  93.           standard 754 for double precision floating-point arithmetic.
  94.  
  95.       ii. On other than VAX or IEEE machines, run the original math 
  96.           library, formerly "/usr/lib/libm.a", now "/usr/lib/libom.a", if
  97.       nothing better is available.
  98.  
  99.      iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s",
  100.           "VAX/tan.s" and "VAX/atan2.s" are different from those in
  101.           "IEEE/trig.c" and "IEEE/atan2.c".  The VAX assembler code uses the
  102.           true value of pi to perform argument reduction, while the C code uses
  103.           a machine value of PI (see "IEEE/trig.c").
  104.  
  105.  
  106. 2. A computer system that conforms to IEEE standard 754 should provide 
  107.                 sqrt(x),
  108.                 drem(x,p), (double precision remainder function)
  109.                 copysign(x,y),
  110.                 finite(x),
  111.                 scalb(x,N),
  112.                 logb(x) and
  113.                 rint(x).
  114.    These functions are either required or recommended by the standard.
  115.    For convenience, a (slow) C implementation of these functions is 
  116.    provided in the file "IEEE/support.c".
  117.  
  118.    Warning: The functions in IEEE/support.c are somewhat machine dependent.
  119.    Some modifications may be necessary to run them on a different machine.
  120.    Currently, if compiled with a suitable flag, "IEEE/support.c" will work
  121.    on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile"
  122.    in this directory). Invoke the C compiler thus:
  123.  
  124.         cc -c -DVAX IEEE/support.c              ... on a VAX, D-format
  125.         cc -c -DNATIONAL IEEE/support.c         ... on a National 32000
  126.         cc -c  IEEE/support.c                   ... on other IEEE machines,
  127.                                                     we hope.
  128.  
  129.    Notes: 
  130.       1. Faster versions of "drem" and "sqrt" for IEEE double precision
  131.          (coded in C but intended for assembly language) are given at the
  132.          end of "IEEE/support.c" but commented out since they require certain
  133.          machine-dependent functions.
  134.  
  135.       2. A fast VAX assembler version of the system supported functions
  136.          copysign(), logb(), scalb(), finite(), and drem() appears in file 
  137.          "VAX/support.s".  A fast VAX assembler version of sqrt() is in
  138.          file "VAX/sqrt.s".
  139.  
  140. 3. Two formats are supported by all the standard elementary functions: 
  141.    the VAX D-format (56-bit precision), and the IEEE double format 
  142.    (53-bit precision).  The cbrt() in "IEEE/cbrt.c" is for IEEE machines 
  143.    only. The functions in files that end with ".s" are for VAX computers 
  144.    only. The functions in files that end with ".c" (except "IEEE/cbrt.c")
  145.    are for VAX and IEEE machines. To use the VAX D-format, compile the code 
  146.    with -DVAX; to use IEEE double format on various IEEE machines, see 
  147.    "Makefile" in this directory). 
  148.  
  149.     Example:
  150.         cc -c -DVAX sin.c               ... for VAX D-format
  151.  
  152.        Warning: The values of floating-point constants used in the code are
  153.                 given in both hexadecimal and decimal.  The hexadecimal values
  154.                 are the intended ones. The decimal values may be used provided 
  155.                 that the compiler converts from decimal to binary accurately
  156.                 enough to produce the hexadecimal values shown. If the
  157.                 conversion is inaccurate, then one must know the exact machine 
  158.                 representation of the constants and alter the assembly
  159.                 language output from the compiler, or play tricks like
  160.                 the following in a C program.
  161.  
  162.                         Example: to store the floating-point constant 
  163.  
  164.                              p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)
  165.  
  166.                         on a VAX in C, we use two longwords to store its 
  167.                         machine value and define p1 to be the double constant 
  168.                         at the location of these two longwords:
  169.  
  170.                         static long  p1x[] = { 0x3abe3d78, 0x066a67e1};
  171.                         #define      p1      (*(double*)p1x)
  172.  
  173.     Note:  On a VAX, some functions have two codes. For example, cabs() has
  174.        one implementation in "IEEE/cabs.c", and another in "VAX/cabs.s". 
  175.            In this case, the assembly language version is preferred.
  176.  
  177.  
  178. 4. Accuracy. 
  179.  
  180.             The errors in expm1(), log1p(), exp(), log(), cabs(), hypot()
  181.             and cbrt() are below 1 ULP (Unit in the Last Place).
  182.  
  183.             The error in pow(x,y) grows with the size of y. Nevertheless,
  184.             for integers x and y, pow(x,y) returns the correct integer value 
  185.             on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that 
  186.             x to the power of y is representable exactly.
  187.  
  188.             cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
  189.             about 3 ULPs. 
  190.  
  191.             For trigonometric and inverse trigonometric functions: 
  192.  
  193.                 Let [trig(x)] denote the value actually computed for trig(x),
  194.  
  195.                 1) Those codes using the machine's value PI (true pi rounded):
  196.                    (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c)
  197.  
  198.                    The errors in [sin(x)], [cos(x)], and [atan(x)] are below 
  199.                    1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and 
  200.                    atan(x)*PI/pi respectively, where PI is the machine's
  201.                    value of pi rounded. [tan(x)] returns tan(x*pi/PI) within
  202.                    about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)] 
  203.                    return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
  204.                    respectively to similar accuracy.
  205.  
  206.  
  207.                 2) Those using true pi (for VAX D-format only):
  208.                    (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and
  209.                    atan.c)
  210.  
  211.                    The errors in [sin(x)], [cos(x)], and [atan(x)] are below
  212.                    1 ULP. [tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)] 
  213.                    have errors below about 2 ULPs. 
  214.  
  215.  
  216.             Here are the results of some test runs to find worst errors on 
  217.             the VAX :
  218.  
  219.     tan   :  2.09 ULPs          ...1,024,000 random arguments (machine PI)
  220.     sin   :  .861 ULPs          ...1,024,000 random arguments (machine PI)
  221.     cos   :  .857 ULPs          ...1,024,000 random arguments (machine PI)
  222.     (compared with tan, sin, cos of (x*pi/PI))
  223.  
  224.     acos  :  2.07 ULPs          .....200,000 random arguments (machine PI)
  225.     asin  :  2.06 ULPs          .....200,000 random arguments (machine PI)
  226.     atan2 :  1.41 ULPs          .....356,000 random arguments (machine PI)
  227.     atan  :  0.86 ULPs          ...1,536,000 random arguments (machine PI)
  228.     (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
  229.  
  230.     tan   :  2.15 ULPs          ...1,024,000 random arguments (true pi)
  231.     sin   :  .814 ULPs          ...1,024,000 random arguments (true pi)
  232.     cos   :  .792 ULPs          ...1,024,000 random arguments (true pi)
  233.     acos  :  2.15 ULPs          ...1,024,000 random arguments (true pi)
  234.     asin  :  1.99 ULPs          ...1,024,000 random arguments (true pi)
  235.     atan2 :  1.48 ULPs          ...1,024,000 random arguments (true pi)
  236.     atan  :  .850 ULPs          ...1,024,000 random arguments (true pi)
  237.  
  238.     acosh :  3.30 ULPs          .....512,000 random arguments
  239.     asinh :  1.58 ULPs          .....512,000 random arguments
  240.     atanh :  1.71 ULPs          .....512,000 random arguments  
  241.     cosh  :  1.23 ULPs          .....768,000 random arguments
  242.     sinh  :  1.93 ULPs          ...1,024,000 random arguments
  243.     tanh  :  2.22 ULPs          ...1,024,000 random arguments
  244.     log10 :  1.74 ULPs          ...1,536,000 random arguments
  245.     pow   :  1.79 ULPs          .....100,000 random arguments, 0 < x, y < 20.
  246.         
  247.     exp   :  .768 ULPs          ...1,156,000 random arguments
  248.     expm1 :  .844 ULPs          ...1,166,000 random arguments
  249.     log1p :  .846 ULPs          ...1,536,000 random arguments
  250.     log   :  .826 ULPs          ...1,536,000 random arguments
  251.     cabs  :  .959 ULPs          .....500,000 random arguments
  252.     cbrt  :  .666 ULPs          ...5,120,000 random arguments
  253.  
  254.  
  255. 5. Speed.
  256.  
  257.         Some functions coded in VAX assembly language (cabs(), hypot() and
  258.     sqrt()) are significantly faster than the corresponding ones in 4.2BSD.
  259.         In general, to improve performance, all functions in "IEEE/support.c"
  260.         should be written in assembly language and, whenever possible, should
  261.     be called via short subroutine calls.
  262.  
  263.  
  264. 6. j0, j1, jn.
  265.  
  266.         The modifications to these routines were only in how an invalid
  267.         floating point operations is signaled.
  268.